home *** CD-ROM | disk | FTP | other *** search
/ Eagles Nest BBS 8 / Eagles_Nest_Mac_Collection_Disc_8.TOAST / Developer Tools⁄Additions / MacScheme20 / Mathlib / chebpoly.scm < prev    next >
Encoding:
Text File  |  1989-04-27  |  7.3 KB  |  209 lines  |  [TEXT/????]
  1. ;;; 12/23/87 (mh) modified CHEB-ECON to return original poly if not trimmed
  2. ;;; $Header: chebpoly.scm,v 1.3 88/01/26 06:58:01 GMT gjs Exp $
  3. ;;; Chebyshev polynomial routines
  4.  
  5. (if-mit
  6.  (declare (usual-integrations = + - * /
  7.                  zero? 1+ -1+
  8.                  ;; truncate round floor ceiling
  9.                  sqrt exp log sin cos)))
  10.  
  11.  
  12. ;;; We define a function that returns a stream of Chebyshev polynomials,
  13. ;;;  using the recurrence relation
  14. ;;;           T[0] = 1, T[1] = x, T[n] = 2xT[n-1] - T[n-2]
  15.  
  16. (define (chebyshev-polynomials)
  17.   (define s
  18.     (cons-stream
  19.       '(1)
  20.       (cons-stream
  21.         '(0 1)
  22.         (map-streams (lambda (p1 p2)
  23.                        (sub-polys (mul-polys '(0 2) p1) p2))
  24.                      (tail s)
  25.                      s))))
  26.   s)
  27.  
  28. ;;; The following procedure returns the Nth Chebyshev polynomial
  29.  
  30. (define (cheb-poly n)
  31.   (let ((s (chebyshev-polynomials)))
  32.     (stream-ref s n)))
  33.  
  34.  
  35. ;;; In the following, we define a CHEB-EXP to be an expansion in
  36. ;;;  Chebyshev polynomials. As with polynomials, the expansion is
  37. ;;;  written as a list of coefficients (a0 a1 ... aN), and represents
  38. ;;;  the formal series a0*T[0](x) + ... + aN*T[N](x).
  39.  
  40. ;;; The following procedure returns a stream of scaled Chebyshev
  41. ;;;  expansions: precisely, it returns the expansions corresponding
  42. ;;;  to the sequence 1, x, 2x^2, 4x^3, ..., 2^(n-1)x^n ... Note that the
  43. ;;;  general form doesn't cover the first term. What is generally wanted
  44. ;;;  is the expansion corresponding to x^n; the power of 2 is thrown in
  45. ;;;  so that the resulting expansion is exact in integer arithmetic.
  46.  
  47. (define (scaled-chebyshev-expansions)
  48.   (define (2x cheb-exp)
  49.     (let ((t1 (cdr cheb-exp))
  50.           (t2 (append (list 0) cheb-exp))
  51.           (t3 (list 0 (car cheb-exp))))
  52.       (add-polys t1 (add-polys t2 t3))))
  53.   (define s
  54.     (cons-stream
  55.       '(1)
  56.       (cons-stream
  57.         '(0 1)
  58.         (map-stream 2x (tail s)))))
  59.   s)
  60.  
  61.  
  62. ;;; For convenience, we also provide the non-scaled Chebyshev expansions
  63.  
  64. (define (chebyshev-expansions)
  65.   (letrec ((s (cons-stream ; s = {1 1 2 4 8 16 ...}
  66.                 1
  67.                 (cons-stream
  68.                   1
  69.                   (map-stream (lambda(x) (+ x x)) (tail s)))))
  70.            (c (scaled-chebyshev-expansions)))
  71.     (map-streams (lambda (factor poly)
  72.                    (scale-poly (/ 1 factor) poly))
  73.                  s
  74.                  c)))
  75.  
  76.  
  77. ;;; Convert from polynomial form to Chebyshev expansion
  78.  
  79. (define (poly->cheb-exp p)
  80.   (let ((n (length p)))
  81.     (let ((cheb (first-n-of-stream (chebyshev-expansions) n)))
  82.       (reduce
  83.         add-polys
  84.         (map scale-poly p cheb)))))
  85.  
  86. ;;; Convert from Chebyshev expansion to polynomial form
  87.  
  88. (define (cheb-exp->poly p)
  89.   (let ((n (length p)))
  90.     (let ((cheb (first-n-of-stream (chebyshev-polynomials) n)))
  91.       (reduce
  92.         add-polys
  93.         (map scale-poly p cheb)))))
  94.  
  95. ;;; Given a Cheb-expansion and an error criterion EPS, trim the tail of
  96. ;;;  those coefficients whose absolute sum is <= EPS.
  97.  
  98. (define (trim-cheb-exp cheb eps)
  99.   (let ((r (reverse cheb)))
  100.     (let loop ((sum (abs (car r))) (r r))
  101.       (if (= (length r) 1)
  102.           (if (<= sum eps) '(0) r)
  103.           (if (> sum eps)
  104.               (reverse r)
  105.               (loop (+ sum (abs (cadr r))) (cdr r)))))))
  106.  
  107.  
  108. ;;; The next procedure performs Chebyshev economization on a polynomial p
  109. ;;;  over a specified interval [a,b]: the returned polynomial is guaranteed
  110. ;;;  to differ from the original by no more than eps over [a, b].
  111.  
  112. (define (cheb-econ p a b eps)
  113.   (let ((q (poly-domain->canonical p a b)))
  114.     (let ((r (poly->cheb-exp q)))
  115.       (let ((s (trim-cheb-exp r eps)))
  116.         (if (= (length s) (length r)) ;nothing got trimmed
  117.             p
  118.             (let ((t (cheb-exp->poly s)))
  119.               (poly-domain->general t a b)))))))
  120.  
  121.  
  122. ;;; Return the root-list of the Nth Chebyshev polynomial
  123.  
  124. (define (cheb-root-list n)
  125.   (define (root i)
  126.     (if (and (odd? n) (= i (/ (- n 1) 2)))
  127.         0
  128.         (- (cos (/ (* (+ i (/ 1 2)) pi) n)))))
  129.   (let loop ((i 0))
  130.     (if (= i n)
  131.         '()
  132.         (cons (root i)
  133.               (loop (+ i 1))))))
  134.  
  135. ;;; This procedure accepts an integer n > 0 and a real x, and returns
  136. ;;;  a list of the values T[0](x) ... T[n-1](x). If an optional third
  137. ;;;  argument is given as the symbol 'HALF, then the first value in the
  138. ;;;  list will have the value 0.5; otherwise it has the value 1.
  139.  
  140. (define (first-n-cheb-values n x . optionals)
  141.   (let ((first (if (and (not (null? optionals))
  142.                         (eq? (car optionals) 'HALF))
  143.                    (/ 1 2)
  144.                    1)))
  145.     (cond ((< n 1) '())
  146.           ((= n 1) (list first))
  147.           ((= n 2) (list first x))
  148.           (else (let loop ((ans (list x first)) (a 1) (b x) (count 2))
  149.                   (if (= count n)
  150.                       (reverse ans)
  151.                       (let ((next (- (* 2 x b) a)))
  152.                         (loop (cons next ans) b next (+ count 1)))))))))
  153.  
  154.  
  155. ;;; The following procedure evaluates a cheb-expansion directly, in a
  156. ;;;  manner analogous to Horner's method for evaluating polynomials --
  157. ;;;  actually, it isn't as analogous as it should be, and should
  158. ;;;  probably be replaced by Clenshaw's method which truly is like
  159. ;;;  Horner's.
  160.  
  161. (define (cheb-exp-value cheb x)
  162.   (let ((n (length cheb)))
  163.     (let ((vals (first-n-cheb-values n x)))
  164.       (reduce + (map * cheb vals)))))
  165.  
  166.  
  167. ;;; This procedure generates a Chebyshev expansion to a given order N,
  168. ;;;  for a function f specified on an interval [a,b]. The interval is
  169. ;;;  mapped onto [-1,1] and the function approximated is g, defined on
  170. ;;;  [-1,1] to behave the same as f on [a,b].
  171. ;;; Note: the returned list of coefficients is N long, and is associated
  172. ;;;  with Chebyshev polynomials from T[0] to T[N-1].
  173.  
  174. (define (generate-cheb-exp f a b n)
  175.   (if (<= b a)
  176.       (error "Bad interval in GENERATE-CHEB-EXP"))
  177.   (let ((interval-map    ;map [-1,1] onto [a,b]
  178.           (let ((c (/ (+ a b) 2))
  179.                 (d (/ (- b a) 2)))
  180.             (lambda (x) (+ c (* d x))))))
  181.     (let ((roots (cheb-root-list n))
  182.           (polys (first-n-of-stream (chebyshev-polynomials) n)))
  183.       (let ((vals (map f (map interval-map roots))))
  184.         (let loop ((coeffs '()) (i 0))
  185.           (if (= i n)
  186.               (reverse coeffs)
  187.               (let ((chebf (lambda(x) (poly-value (list-ref polys i) x))))
  188.                 (let ((chebvals (map chebf roots)))
  189.                   (let ((sum (reduce + (map * vals chebvals))))
  190.                     (let ((term (if (zero? i) (/ sum n) (/ (* 2 sum) n))))
  191.                       (loop (cons term coeffs) (+ i 1))))))))))))
  192.  
  193.  
  194. ;;; This procedure accepts a function f, an interval [a,b], a number
  195. ;;;  N of points at which to base a Chebyshev interpolation, and an
  196. ;;;  optional precision EPS. The method is to map f onto the canonical
  197. ;;;  interval [-1,1], generate the Chebyshev expansion based on the
  198. ;;;  first N Chebyshev polynomials interpolating at the roots of the
  199. ;;;  N+1st Chebyshev polynomial T[N]. If EPS has been specified, an
  200. ;;;  economization is performed at this point; otherwise not. Finally,
  201. ;;;  the Chebyshev expansion is reconverted to a polynomial and mapped
  202. ;;;  back onto the original interval [a,b].
  203.  
  204. (define (generate-approx-poly f a b n . optionals)
  205.   (let ((eps (if (null? optionals) false (car optionals))))
  206.     (let ((p (generate-cheb-exp f a b n)))
  207.       (let ((pp (if eps (trim-cheb-exp p eps) p)))
  208.         (poly-domain->general (cheb-exp->poly pp) a b)))))
  209.